rm(list = ls(all = TRUE))
library(foreign)
library(mediation)

d <- read.dta("KS_replication_european_monarchies.dta")

## Causal Mediation
med.dat <- function(dat, outcome){
    ids <- sort(unique(dat$countrycode))
    ndat <- c()
    for(i in 1:length(ids)){
        temp <- dat[dat$countrycode == ids[i], ]
        temp$lag_tenure <- NA
        temp$lag_deposed <- NA
        for(j in 2:nrow(temp)){
            temp$lag_tenure[j] <- temp$tenure[j-1]
            if(outcome == "ks") temp$lag_deposed[j] <- temp$our_deposed[j-1]
            if(outcome == "dow") temp$lag_deposed[j] <- temp$dow_deposed[j-1]
        }
        ndat <- rbind(ndat, temp)
    }
    dsub <- ndat[ndat$lag_deposed == 0, ]
    qt <- quantile(dsub$lag_tenure, c(0.25, 0.5, 0.75), na.rm = TRUE)
    dsub$tenure_qt <- ifelse(dsub$lag_tenure < qt[1], 1, ifelse(dsub$lag_tenure >= qt[1] & dsub$lag_tenure < qt[2], 2, ifelse(dsub$lag_tenure  >= qt[2] & dsub$lag_tenure < qt[3], 3, 4))) 
    if(outcome == "dow") dsub <- dsub[is.na(dsub$dow_deposed) == FALSE, ] ## Use when analyzing dow_deposed as the outcome
    return(dsub)
}

out.med <- function(med.fit, out.fit){

    med.out12 <- mediate(med.fit, out.fit, treat = "tenure_qt", mediator = "primogeniture", control.value = 1, treat.value = 2, sims = 1000)
    med.out13 <- mediate(med.fit, out.fit, treat = "tenure_qt", mediator = "primogeniture", control.value = 1, treat.value = 3, sims = 1000)
    med.out14 <- mediate(med.fit, out.fit, treat = "tenure_qt", mediator = "primogeniture", control.value = 1, treat.value = 4, sims = 1000)
    med.out23 <- mediate(med.fit, out.fit, treat = "tenure_qt", mediator = "primogeniture", control.value = 2, treat.value = 3, sims = 1000)
    med.out24 <- mediate(med.fit, out.fit, treat = "tenure_qt", mediator = "primogeniture", control.value = 2, treat.value = 4, sims = 1000)
    med.out34 <- mediate(med.fit, out.fit, treat = "tenure_qt", mediator = "primogeniture", control.value = 3, treat.value = 4, sims = 1000)


    outMULT <- data.frame(rbind(c(med.out12$z.avg.ci[1], med.out12$z.avg, med.out12$z.avg.ci[2]),
                                c(med.out12$d.avg.ci[1], med.out12$d.avg, med.out12$d.avg.ci[2]),
                                c(med.out13$z.avg.ci[1], med.out13$z.avg, med.out13$z.avg.ci[2]),
                                c(med.out13$d.avg.ci[1], med.out13$d.avg, med.out13$d.avg.ci[2]),
                                c(med.out14$z.avg.ci[1], med.out14$z.avg, med.out14$z.avg.ci[2]),
                                c(med.out14$d.avg.ci[1], med.out14$d.avg, med.out14$d.avg.ci[2]),
                                c(med.out23$z.avg.ci[1], med.out23$z.avg, med.out23$z.avg.ci[2]),
                                c(med.out23$d.avg.ci[1], med.out23$d.avg, med.out23$d.avg.ci[2]),
                                c(med.out24$z.avg.ci[1], med.out24$z.avg, med.out24$z.avg.ci[2]),
                                c(med.out24$d.avg.ci[1], med.out24$d.avg, med.out24$d.avg.ci[2]),
                                c(med.out34$z.avg.ci[1], med.out34$z.avg, med.out34$z.avg.ci[2]),
                                c(med.out34$d.avg.ci[1], med.out34$d.avg, med.out34$d.avg.ci[2])))

    colnames(outMULT) <- c("ci.lower", "pe", "ci.upper")
    outMULT$type <- NA
    outMULT$type[c(1, 3, 5, 7, 9, 11)] <- "ade"
    outMULT$type[c(2, 4, 6, 8, 10, 12)] <- "acme"
    outMULT$treat[1:2] <- 1
    outMULT$treat[3:4] <- 2
    outMULT$treat[5:6] <- 3
    outMULT$treat[7:8] <- 4
    outMULT$treat[9:10] <- 5
    outMULT$treat[11:12] <- 6

    outMULT$pos[outMULT$treat == 1] <- c(0.75, 1.25)
    outMULT$pos[outMULT$treat == 2] <- c(2.75, 3.25)
    outMULT$pos[outMULT$treat == 3] <- c(4.75, 5.25)
    outMULT$pos[outMULT$treat == 4] <- c(6.75, 7.25)
    outMULT$pos[outMULT$treat == 5] <- c(8.75, 9.25)
    outMULT$pos[outMULT$treat == 6] <- c(10.75, 11.25)
    return(outMULT)
}


## Morby Outcome
ddow <- med.dat(d, "dow")
med.fit <- glm(primogeniture ~ tenure_qt, family = binomial("logit"), data = ddow)
out.fit <- glm(dow_deposed ~ primogeniture + tenure_qt, data = ddow, family = binomial("logit"))

eddow <- out.med(med.fit, out.fit)

pdf("figure6left.pdf")
par(mar = c(4, 4.5, 2.5, 0))
plot(eddow$pos, eddow$pe, ylim = c(-0.25, 0.05), type = "n", yaxt = "n", main = "Morby", xaxt = "n", xlab = "Leader Duration (in Years) Quartile Condition Comparison",
     ylab = "", cex.lab = 1.35, cex = 1.25, cex.main = 1.5, frame = FALSE)
points(eddow$pos[eddow$type == "acme"], eddow$pe[eddow$type == "acme"], pch = 18)
points(eddow$pos[eddow$type == "ade"], eddow$pe[eddow$type == "ade"], pch = 20)
abline(h = 0, lty = "dashed")
axis(2, seq(-0.25, 0.05, by = 0.05), seq(-0.25, 0.05, by = 0.05), cex.axis = 1, las = 2)
axis(1, c(1, 3, 5, 7, 9, 11), c("2 vs. 1", "3 vs. 1", "4 vs. 1", "3 vs. 2", "4 vs. 2", "4 vs. 3"), cex.axis = 1.25, tick = FALSE)
for(i in 1:6) lines(c(eddow$pos[eddow$type == "ade"][i], eddow$pos[eddow$type == "ade"][i]), c(eddow$ci.lower[eddow$type == "ade"][i], eddow$ci.upper[eddow$type == "ade"][i]), lty = "dashed")
for(i in 1:6) lines(c(eddow$pos[eddow$type == "acme"][i], eddow$pos[eddow$type == "acme"][i]), c(eddow$ci.lower[eddow$type == "acme"][i], eddow$ci.upper[eddow$type == "acme"][i]))
legend("bottomright", c("ACME", "ADE"), lty = c("solid", "dashed"), cex = 1.15, bty = "n")
mtext("Change in Probability of Deposal", 2, outer = TRUE, line = -1, cex = 1.25)
dev.off()
graphics.off()

## KS Outcome
dks <- med.dat(d, "ks")
med.fit <- glm(primogeniture ~ tenure_qt, family = binomial("logit"), data = dks)
out.fit <- glm(our_deposed ~ primogeniture + tenure_qt, data = dks, family = binomial("logit"))

edks <- out.med(med.fit, out.fit)

pdf("figure6right.pdf")
par(mar = c(4, 3.2, 2.5, 0))
plot(edks$pos, edks$pe, ylim = c(-0.25, 0.05), type = "n", yaxt = "n", main = "Kokkonen and Sundell", xaxt = "n", xlab = "Leader Duration (in Years) Quartile Condition Comparison",
     ylab = "", cex.lab = 1.35, cex = 1.25, cex.main = 1.5, frame = FALSE)
points(edks$pos[edks$type == "acme"], edks$pe[edks$type == "acme"], pch = 18)
points(edks$pos[edks$type == "ade"], edks$pe[edks$type == "ade"], pch = 20)
abline(h = 0, lty = "dashed")
axis(2, seq(-0.25, 0.05, by = 0.05), seq(-0.25, 0.05, by = 0.05), cex.axis = 1, las = 2)
axis(1, c(1, 3, 5, 7, 9, 11), c("2 vs. 1", "3 vs. 1", "4 vs. 1", "3 vs. 2", "4 vs. 2", "4 vs. 3"), cex.axis = 1.25, tick = FALSE)
for(i in 1:6) lines(c(edks$pos[edks$type == "ade"][i], edks$pos[edks$type == "ade"][i]), c(edks$ci.lower[edks$type == "ade"][i], edks$ci.upper[edks$type == "ade"][i]), lty = "dashed")
for(i in 1:6) lines(c(edks$pos[edks$type == "acme"][i], edks$pos[edks$type == "acme"][i]), c(edks$ci.lower[edks$type == "acme"][i], edks$ci.upper[edks$type == "acme"][i]))
legend("bottomright", c("ACME", "ADE"), lty = c("solid", "dashed"), cex = 1.15, bty = "n")
dev.off()
graphics.off()
